arXiv:1502.05802vl [physics.class-ph] 20 Feb 2015 


Inter-sensor propagation delay estimation using sources of opportunity 

Remy Vincent®’*^, Mikael Carmona*^, Olivier Michel*^, Jean-Louis Lacoume*^’’^ 


^CEA, Leti, 38054 Grenoble, Franee 
^ GIPSA-lab, BP 46 F-38402 Grenoble Gedex, Franee 


Abstract 

Propagation delays are intensively used for Structural Health Monitoring or Sensor Network Localization. In this paper, 
we study the performances of acoustic propagation delay estimation between two sensors, using sources of opportunity 
only. Such sources are defined as being uncontrolled by the user (activation time, location, spectral content in time 
and space), thus preventing the direct estimation with classical active approaches, such as TDOA, RSSI and AOA. 
Observation models are extended from the literature to account for the spectral characteristics of the sources in this 
passive context and we show how time-filtered sources of opportunity impact the retrieval of the propagation delay 
between two sensors. A geometrical analogy is then proposed that leads to a lower bound on the variance of the 
propagation delay estimation that accounts for both the temporal and the spatial properties of the sources field. 


1. Introduction to passive estimation 

Propagation delay estimation is a fundamental issue for 
several applications such as Structural Health Monitoring 
(SHM) and Sensor Network Localization (SNL). For in¬ 
stance, structural monitoring can be achieved by relating 
structural deformations to the variations of a propagation 
delay, between two sensors of known location. SNL can 
be achieved from a set of inter-sensor propagation delays 
by using dedicated algorithms (see e.g. EDI) and benefits 
SHM for at least two reasons; firstly it allows to associate 
the measurements (e.g. temperature) to the precise sen¬ 
sor locations, which can be crucial for diagnosis; secondly, 
it enables the retrieval of physical or geometrical param¬ 
eters from the geometry of the sensor network itself {e.g. 
velocity maps in seismology m)- 

Recently, SNL in aerial acoustics was used for an SHM 
application, namely to retrieve the bending of a beam, see 
m- Microphones were distributed at the surface of the 
structure, and identification of the shape of the beam was 
conducted by processing acoustic pressure fields propagat¬ 
ing around the beam, unlike mechanical displacements. 

Furthermore, as in [23], SNL was performed using 
sources of opportunity only. These latter are defined as 
sources that are uncontrolled (unknown positions, acti¬ 
vation and spectral content), thus preventing the direct 
estimation of the inter-sensor distance with classical ac¬ 
tive approaches, such as TDOA, RSSI and AOA, see e.g. 
[3l[llll0llll[2]. Fig. [^displays an example of source of 
opportunity in aerial acoustic: a handclap. In this case, 
the tail of the recorded responses are spatially correlated, 
due to the propagation in the medium. In sectionj^ a tech¬ 
nique is presented to retrieve the Green’s function between 
the two sensing locations from those responses. Exploiting 
sources of opportunity to estimate parameters of a prop¬ 


agation medium is usually referred to as passive identifi¬ 
cation [6], output-only identification m or seismic inter¬ 
ferometry [lisa- Unlike active methods, passive methods 
do not require any additional equipment to the sensor net¬ 
work, nor does it require the structure to be temporarily 
put out-of-order. 



Figure 1: (a) Schematic of the propagation medium. 

Sources si, s2 and s3 are sources of opportu¬ 
nity. Their location and spectral content are 
uncontrolled, (b) Response of a source of op¬ 
portunity (a handclap). 
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Pioneering work in passive identification from Nyquist 
[18] is related to the fluctuation dissipation theorem El 
and establishes the relationship between the two-point cor¬ 
relation of the field and the Green’s function between those 
two locations, provided the source field is white in time 
and space. Major contributions can be found in seismic 
interferometry for the passive imaging of the Earth inte¬ 
rior [H |5l [U |27l [ 12 ] [ 21 ] , acoustics (ultrasounds) [T71I23] 

and electromagnetism [32] for the local imaging of hetero¬ 
geneities. Early contributions to SHM can be found under 
the name output-only structural analysis, with El illus¬ 
trates the passive monitoring of a bridge from its ambient 
mechanical vibrations. More recent work can be found 
alongside with formalization efforts, see e.g. [anziiin]. 

Interestingly, the passive identification protocols have 
proven to be robust and repeatable in various media and 
at various wavelengths. However, the propagation veloci¬ 
ties of mechanical displacement fields prevents to localize 
at the scale of the structure, because of the current tech¬ 
nological limitations. In this regard, acoustic wave propa¬ 
gation was shown in [30] to lead to a millimetric precision 
for the inter-sensors distance estimation. Thus, acoustic 
propagation (in the bandwidth [0.1 — 20] kHz) represents 
a complementary approach for passive structural analysis. 
Nonetheless, compared to the low frequency and narrow 
band signals from long range propagation considered in 
seismic applications, pressure fields in aerial acoustics are 
rather wide-band and high frequency. As a consequence, 
there is a need to evaluate the performances of passive 
estimators in aerial acoustics and to take the evolutions 
of dissipation effects over the considered frequency range 
into account. To the best of our knowledge, such perfor¬ 
mances study was never derived in a systematic way in 
the literature. The closest investigations to our purpose 
were those of [28], in which the impact of instrumentation 
time shifts on travel-time retrieval is studied. The effects 
of the spatial distribution of sources on the retrieval of the 
Green’s function for a free acoustic propagation medium 
were observed and formalized in 0 [221 EH, but did not 
lead to performances assessment. 

In this paper, we propose an estimation theoretic ap¬ 
proach for passive distances retrieval in the context of ho¬ 
mogeneous acoustic propagation. The first contribution is 
the derivation of a passive estimator of inter-sensor dis¬ 
tances, or delays of propagation, that generalizes existing 
equations obtained in an ideal framework (when the op¬ 
portunity sources are white in time and space), see e.g. 
[331 E]. Secondly, a geometrical analogy of the passive 
estimation process is proposed which leads to analytical 
developments. Erom the proposed observation model, a 
Cramer-Rao lower bound is derived for the propagation 
delay estimation using sources of opportunity only. Simu¬ 
lations and real-data experiments illustrate the relevance 
of our results and their applications in a practical context. 

This paper is organized as follows: 

• The passive context is first presented and the iden¬ 


tity that provides an estimator of Green’s functions is 
extended for fields that are not white in time; 

• Three passive estimators are then investigated, for the 
Green’s correlation function, for the Green’s function 
and for the propagation delay. They are studied in 
terms of the operations required to compute them and 
the parameters that influence their performances; 

• A novel observation model is then derived that allows 
to account for both the spatial and temporal proper¬ 
ties of the sources, in the passive estimation process; 

• A lower bound on the mean square error for the pas¬ 
sive propagation delay estimation is presented. 


2. Passive Green’s function retrieval 

Based on the prior work of [ 8 ], three passive estimators 
are presented that account for non whiteness (in time) 
of the source over the considered bandwidth; estimators 
for the Green’s correlation, the Green’s function and the 
propagation delay. 


2.1. Propagation equation and Greenes function 

Let A be a linear acoustic wave propagation medium, 
also assumed isotropic and homogeneous. The propaga¬ 
tion equation that describes fields in fluids and in linear 
media with damping was derived by Stokes, see e.g. [2ai3]. 
The pressure field / is related to the source field s by 


"A 


+ ,A-+„“A 


/(t,x) = s(t,x). 


( 1 ) 


where v is the wave velocity and 77 is the kinematic viscos¬ 
ity, related to damping. Indeed, the above equation dis¬ 
plays a damping term 77 A^/(t,x), which indicates that 
damping increases proportionally to the cube of the fre¬ 
quency. As short to medium range propagation distances 
are considered — between source and sensors, as well as 
between two sensors — high frequency components of the 
source spectra are significant and cannot be overlooked. 

When the propagation operator in Equation 0 is in¬ 
vertible, its kernel is the Green’s function ^(t — to, x, y). It 
is the impulse response in A between locations x and y, as 
the source transmits at time to- The Green’s function al¬ 
lows to relate sources and responses by a filtering process 
generalized to time and space variables. Eor unbounded 
media, it reads 

/(t,x) = // g{t — t'— \x)s{t'^\x)dt'dvL. ( 2 ) 

J 

In the following, the convolution operator * is used 
to describe filtering operations. Subscripts indicate the 
dimension (time/space) along which the convolution is 
performed, such that Equation Q is then equivalent to 
/(t,x) = [g * s](t,x). 

t,U 
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2.2. Green’s correlation and Ward identity 


Pioneering work in [33] elaborates on the fluctuation 
dissipation theorem and shows in acoustics that there ex¬ 
ists a relationship between the two-point correlation of the 
field and the Green’s function between those two locations, 
provided the source field is white in time and space. This 
identity is here extended to sources in aerial acoustics (in 
the audio range [0.1 — 20] kHz) that are not white in time. 
Firstly, the correlation of the field / recorded at x and y 
is defined as 


'lf{u,v,x,y) = E[f{u,x)f*{v,y)], (3) 

where the expectation operator E applies on both the time 
and the space variable. 

In the literature, see for instance [33l nil El [ig, the 
source field is assumed to be stationary and white in both 
time and spatial domains. Its correlation function is 

E[s(t,x)s*(t+ M,y)] =(5(M)(5(x-y). (4) 

In this case, the correlation of the fields recorded is 
called the Green’s correlation function. It is noted 7 ^ and 
is obtained combining Equations ( 11 ,® and 0 , 

Igiu, X, y) = [g * g~] {u, x - y). (5) 

where function g~ is the time-reversed version of e.g. 

In |33l 1121 1?^. an identity is derived from Equation 0 
that relates the Green’s correlation 7 ^ function to the odd 
part of the Green’s function oddf^f] = ^{g — g~). The 
resulting relationship is called Ward identity and accounts 
for constant damping, which means the damping term in 
Equation 0 is proportional to ^/(t,x), such that 

^ 7 ^(^,x-y) = -ioddf^f] (^,x,y). 

where it appears that the damping (7 7 ^ 0 ) is necessary 
to establish the Ward identity. A Ward identity was more 
specifically derived for weakly viscous acoustic media in 
[ 6 ]. Weak viscosity allows to assume g/vX where A is 
the wavelength. Under this assumption, the Ward identity 
can be approximated by 


2.3. Generalization to non temporally white sources 

Unlike in Equation ([^, assume instead that the (spa¬ 
tially extended) source s is white in space but is stationary 
in time. Its time auto-correlation function is 7 s (t) and its 
total correlation function is given by 

E[s(f,x)s*(f + u,y)] =7s(u)(5(x-y). (8) 

Combining Equations 0,0 and ^ results in an ex¬ 
tended interferences formula 


^[f{t,^)f*{t + u,y)] = [ys*g *^g ](u,x-y). (9) 


where the correlation of the source field acts as a filter 
on the Green’s correlation from Equation Combining 
Equations ^ and 0 to the Green’s function estimator 
g{u^ X, y) in Equation 0 results in an estimator that ac¬ 
counts for the spectral properties of the source field. 


g{u,x,y) « todd7 * 5 ](u,x,y). 
g ^ t 


( 10 ) 


As can be seen in Equation (10), sources of opportunity 


allow to retrieve Green’s functions filtered by the correla¬ 
tion of the source, which is related to the spectral content 
of the source itself. Eor a source that is white in time and 
space, the Green’s function could theoretically be retrieved 
with an arbitrarily small estimation residue, depending on 
the observation duration and the sensors noise only. 


3. Passive parameter estimation 
3.1. Observation model 

Physical parameters characterizing the propagation 
medium as presented in the introduction may be inferred 
from the Green’s function. Here, highlight is made on the 
passive retrieval of inter-sensor propagation delays. The 
minimum propagation delay is defined as the time required 
for a wave to travel from a position x to a position y by the 
shortest path. If the velocity of the wave in the medium 
is known, an inter-sensor distance can be easily retrieved 
from this minimum propagation delay. 

Assume that the dimensions of the medium X are much 
larger than the wavelength; then quasi-optical ray prop¬ 
agation is considered [ig. The Green’s function can be 
modeled as. 


^3 

dt"^ 


lg{u,x-y) 


1 

7] 


Odd [5] (u,x,y). 


CX3 

( 6 ) 9 {t,x,y) = '^aiS{t-ti), with |ai| < |aj|, i < j, ( 11 ) 


In other words, the Ward identity provides an estimator 
^(i4, X, y), for the Green’s function between locations x 
and y. Because the Green’s function is causal, retrieving 
its odd part suffices to retrieve it entirely, 

g3 

g{u,X,y) =-^jg{u,x-y), for u>0. (7) 


where the amplitudes decay due to damping and geomet¬ 
rical attenuation. In the following, detection/estimation 
performances are presented for the passive minimum prop¬ 
agation delay estimation, ti. Erom model ( 0 , a simple 
estimator U for ti can be 

= argmax { 5 ((t,x,y)}. ( 12 ) 
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In the passive context, it was shown in Equation (10) mum propagation delay in Equation (12) is extracted from 


that the Green’s function is estimated as a function of the 
source auto-correlation function. If the Green’s function 
is retrieved at short times, it emulates the response of the 
same medium but with no bounds, thus involving the first 
propagation delay only. In this case, an observation model 
for the retrieved Green’s function can be 


git) = - ti) + wit), 


(13) 


the Green’s function estimated using sources of opportu¬ 
nity. This estimation process is thus once again affected 
by the spectral content of the source field, as expressed in 


Equations (10). In the following, the content of the time 


frequency spectrum and the spatial frequency spectrum 
of the source are in turn discussed, with respect to their 
impact on passive estimation performances. 


where w is an additive estimation noise, assumed to be 
white Gaussian, with variance for the sake of argu¬ 
mentation. js is the normalized correlation of the source. 


3.2. Theoretical performances 
3.2.1. Detection 

In the active context, propagation delay estimation can 
be performed by locating the maximum amplitude of 
the inter-correlation between the reference source signal 
and the propagated signal. Interestingly, the observation 
model (13) is very similar to the one encountered in an 
active context. In the passive context however, the source 
is unknown but is expected to be as white as possible 
^ ( 5 ), to lead to optimal performances. Eor this rea¬ 
son, the propagation delay estimator (El is thus justified 
in the passive context as well. 

Erom model (13), taking the maximum amplitude of g{t) 
results in ^4, unless the amplitude of the estimation noise 
gets higher than A at some point. The probability of this 
happening is a standard probability detection error, illus¬ 
trated in Eig. see m for an example. The probability 
that the amplitude of the noise (assumed Gaussian) gets 
higher than A writes 


3.2.2. Estimation 

Let the source be white over the bandwidth B. The pas¬ 
sive retrieval of the Green’s function at short times (in¬ 
volving the first propagation delay only) is thus derived 
from model (El as 

g{t) = A sine(^B{t — ti))w{t) (15) 

where the estimation residue w{t) is modeled Gaussian, 
with power spectral density constant equal to N^j on the 
bandwidth B. Let the sampling rate be A = l/{2Bs). 
Interestingly, this model applies for both the active and the 
passive context, tough with different levels of the signal of 
interest and of the noise. Eor this model, the lower bounds 
on the variance for the propagation delay estimation can 
be found in the literature, for instance in m, such that 


E[(ii - ti)^] > 


N„ 




(16) 


where is shown to be a SNR in (31] and ^tt^BBs 
is identified as the mean square bandwidth of the signal, 
such as in m for instance. 


P[«;[n]>A] = hl-erf(Afi)), (14) 

z Z a 

where erf(x) = 2/yEr exp(—t^)dt is the error function. 
The probability increases as the ratio A/a decreases. Be¬ 
cause the noise samples are assumed i.i.d., the probability 
of misdetection is proportional to the length of the signal. 



Figure 2: Model for Green’s function g(t) = AS(t — ti) + 
w(t) and probability error of extracting the 
pulse location. 

Recall that in the passive context, the Green’s function 
is estimated using sources of opportunity. Thus, the mini- 


3.3. Illustration for time-filtered signals 

In this section, the previous theoretical performances are 
illustrated on a real-data experiment, which consists of the 
passive estimation of the acoustic propagation delay be¬ 
tween two microphones in a bounded environment. Low- 
pass filters with various cut-off frequencies were applied to 
the recorded signals and the estimation performances were 
computed for each of these band-limited signals. 

Sources of opportunity can be pulses emitted at uncon¬ 
trolled times, locations and with uncontrolled energies. In 
this case, it was shown in [H |33] that the tail of those 
impulse responses, called coda waves [T], are the part of 
the signals that can be used to estimate the Green’s func¬ 
tion between the two sensing locations. Here, handclaps 
randomly distributed in the room, were used as sources of 
opportunity. In the experiment, two omni-directional mi¬ 
crophones were placed 48 cm away from each other, in a 
room with dimensions 10 mx 6 mx3.5m. Sensors recorded 
time series x[n] and y[n] at sampling frequency 44.1 kHz, 
see Eigure[^ As in [^, the mixing time was used to set 
the beginning of the coda waves. A threshold was used to 
define the end of the coda waves, as the SNR decays in 
time due to damping and geometrical attenuation. 


4 













The Green’s correlation was estimated from inter¬ 
correlations of coda waves, using the protocol described in 
Appendix I Appendix A[ Extracted codas were lowpass fil¬ 
tered with various cut-off frequencies. For each bandwidth, 
250 inter-correlations of 5-ms segments of the coda waves 
were averaged. Segments overlapped with each others 
and were smoothed by a Blackman apodization window. 
The Green’s functions were estimated from the retrieved 
Green’s correlations, by the application of the Ward iden¬ 
tity ( prl ). Results are displayed in Fig. in which the 
reference (ref) indicates the expected position of the prop¬ 
agation delay. 



Figure 3: (left) Green’s correlation estimated from noise 
correlation. Source fields were lowpass filtered, 
with cut-off frequencies 22 kHz (a), 18 kHz (b), 
9 kHz (c), 5 kHz (d) and 2 kHz(e). (right) Es¬ 
timated Green’s functions. The (ref) function 
locates the propagation delay at 1.40 ms. No 
equivalent reference for the Green’s correlation 
is available in this experiment. 


In this experiment, handclaps were found to have en¬ 
ergy up to 15 kHz. When lowpass filtering the signals, the 
SNR is thus expected to increase with the bandwidth until 
this threshold of 15 kHz is reached, after which only mea¬ 
surement noise remains. The performances for the passive 
propagation delay estimation are thus expected to improve 
up until the cut-off frequency reaches 15 kHz. 

In region (a), increasing the bandwidth leads to a de¬ 
creasing number of misdetections, which in turn results 
in a smaller mean standard deviation. This result was 
predictable from Equation (14), that quantifies the proba¬ 
bility of misdetection as a monotonic function of the SNR. 
As expected, the SNR increases with the bandwidth. 

In region (b) (the cut-off frequency is in the range [7—15] 
kHz), the probability of misdetection as defined in ( p^ de¬ 
creases below 0.001. The lower bound that was derived in 
Equation (16) is overruled by the technical limitations of 
the experimental setup. More precisely, two factors come 
into play: the resolution of the estimation, due to the fact 
that the estimated propagation delays are multiples of the 



Figure 4: Mean standard deviation of the estimated 
propagation delay ti. (dotted line) Resolu¬ 
tion due to the sampling frequency and the 
finite number of iterations, at 9.7 x 10“^ ms. 
(thick line) Performances computed from all 
propagation delays that were estimated, (fine 
line) Performances computed excluding misde¬ 
tections. Discrepancies in the hatched zones 
are clear consequences of misdetection. 


sampling period T^, and the fact that the statistics are 
estimated over a finite number of iterations (1000 in this 
experiment). Fig. illustrates that mostly two neigh¬ 
boring locations are estimated with high probabilities. In 
other words, estimating any other locations than the sam¬ 
ple n or n -|- 1 occurs with a probability that is set by a 
given threshold e. Then the probability that the propaga¬ 
tion delay is contained inside a segment of length 3Ts is 
greater than l-2e. Furthermore, if the probability distribu¬ 
tion function of fho propagation delay is modeled by a 
normal distribution with mean ti and standard deviation 
cr, the following relationship holds: erf () > 1 — 2e. 

In the experimental setup, = 1/44100 and e = 10“^. 
The standard deviation thus meets a resolution bound at 
9.7 X 10“^ ms, which is validated by the experimental re¬ 
sults. For comparison, application of Equation (16) using 
the data in Fig. (c) to estimate the estimation noise 
level over the bandwidth leads to a theoretical bound for 
the standard deviation of 10“^ ms at 9 kHz. 



Figure 5: Distribution of the estimated propagation de¬ 
lay in aerial acoustics. Two samples are esti¬ 
mated with strong probabilities, resulting in a 
relationship between the variance of the esti¬ 
mator and the sampling rate. 
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In region (c) (the bandwidth is larger than 15 kHz), mis- 
detections become more likely, leading to an increase of the 
mean standard deviation. This is interpreted as a severe 
decrease of the SNR for bandwidths that are larger than 
15 kHz, since handclaps have no significant energy above 
this frequency and only measurement noise is recorded. 

As a summary, detecting the propagation delay is more 
likely to fail when the cut-off frequency is outside of the 
band [7 — 15] kHz, which directly results in a strong in¬ 
crease of the standard deviation of the propagation delay 
estimation. On the opposite, when the cut-off frequency is 
inside of the band [7— 15] kHz, the experimental standard 
deviation is seemingly bounded, although it is only due 
of the finite sampling rate and the finite number of itera¬ 
tions of the estimation. Nevertheless, in the case of aerial 
acoustics, this limitation allows to retrieve inter-sensor dis¬ 
tances with a standard deviation that is comparable to the 
dimensions of the microphones 1 cm). 


plane wave sources with a given distribution of incidence 
angles poi see Fig. 



l-L : {A = cste} 



4. An interpretation of the spatial distribution of 
the source field 

In the previous experiment, coda waves were used to 
perform inter-sensor propagation delay estimation in a 
passive context. In this section, a model for coda waves 
is briefly presented via a geometrical analogy, in order to 
easily account for the spatial distribution of the source 
field. Dissipation is not accounted for. We show how the 
spatial distribution impacts the passive propagation de¬ 
lay estimation. The study is illustrated with simulations 
staging space-filtered sources. 

4 . 1 . Towards a CRLB in the passive context 

In this analogy, the propagation medium is assumed 
boundary free. We show that any field that induces time 
differences of arrival at the sensor set could be replaced by 
a set of plane waves with particular incidence angles, as 
suggested in [24]. This allows to model sources of oppor¬ 
tunity (such as the coda wave) with a strong highlight on 
their spatial characteristics. 

Consider a point-wise source at location s in the 
medium, with spherical wavefront. The sensors receive 
the field at different times, delayed by the time difference 
of arrival A that is due to the source-sensors geometry. It 
is shown in [Appendix B| that for any point-wise source, 
there exists an injective transformation T that maps the 
source location s to an incidence angle 6 >eg, so that a plane 
wave impinging on the sensor set given this incidence angle 
would preserve the time difference of arrival, 

T:S^0eq,S.i. A{s)=A{0eq). (17) 


Figure 6: Boundary-free medium analogy. For any time 
difference of arrival A induced by a point-wise 
source, there exists an incidence angle Oeq such 
that a plane wave impinging on the sensor set 
preserves this delay. 


To go further, consider a set of sources of opportunity 
that are i.i.d. in time and space. All fields thus have the 
same time correlation function 7 s (t) and impinge on the 
set of sensors with incidences angles 0 that have the same 
distribution pQ. Denote /(t,x|^^) the field measured at 
position X, due to the source field impinging with the inci¬ 
dence angle Oi and amplitude 5 (t, Oi). The inter-correlation 
of the fields measured at x and y, depends on the associ¬ 
ated time difference of arrival A^, illustrated in Fig. In 
this analogy, the choice is made not to account for attenu¬ 
ation and to focus only on the distribution of the sources, 
in which case the correlation of the fields is 


E[/(i,x|6>j)/*(i + M,y|6>j)] = 7 s(m-A i)(5(6>j - 6>j). (18) 


where 7 ^ is normalized correlation of the source. 

When using sources of opportunity, the field amplitudes 
of interest, e.g. coda waves, are random processes whose 
statistics allow to estimate the Green’s correlation of the 
medium. Everything else is treated as a perturbation, e.g. 


a measurement noise. Therefore, property (18) is funda¬ 


mental to derive an observation model for the Green’s cor¬ 
relation function. Gombined with the correlator in |Ap-| 
jpendix A| the Green’s correlation is estimated by 


p27T 

79 W = / P 0 {O)'ys{t - ticos{0))dO + e{t). (19) 

Jo 


Combining this geometrical property with the Huygens- 
Fresnel principle, any random source field can be modeled 
by an equivalent set of plane wave sources. Furthermore, 
in this analogy, a coda wave is thus replaced by a set of 


where 7 s (t) is the auto-correlation function of the source. 
It is estimated from signals that contain measurements 
noises and that have finite-length, see [Appendix A| Corre¬ 
sponding estimation errors are represented by the additive 


6 













noise e{t). For the sake of simplicity, e{t) is approximated 
centered Gaussian, see [31] for some justifications. 


Equation (19) provides a novel observation model, di¬ 


rectly on the Green’s correlation function. The Fisher in¬ 
formation on the minimum propagation delay ti writes 


tM = 4 / [ 

^ JR^JO 


I cos(^)^^(t—ti cos{0))d0 


dt. 


As in m, a Cramer-Rao Lower Bound for the passive 
context can thereby be derived as 




( 20 ) 


In order to conduct the calculus in Equation (10), the 


strong assumption of spatial whiteness of the source field 
was assumed. However, spatial properties of the sources 
are expected to impact the performances of passive estima¬ 
tors, just as temporal properties do. In the following, the 
new observation model (19) is used to run simulations and 


the passive Cramer-Rao Lower Bound (20) is computed 
and compared to the obtained results. 


4.2. Illustrations with space-filtered sources 

From the model introduced in section |4.H four setups 
are investigated in order to study typical and realistic 
pathologies for the spatial distribution of the source field. 

(1) The first setup involves uniformity of the source field; 

(2) the second setup stages setup 1 with a small count of 
sources; (3) the third setup stages setup 1 with reduced 
spatial bandwidth (solid angle containing the impinging 
wave vectors). (4) Finally, the fourth setup investigates 
what happens when the source field is not uniformly dis¬ 
tributed. Eventually, performances are compared to the 
theoretical bounds obtained in this paper. 


4.2.1. Simulations 

Time differences of arrival are confined to [—1,1] s, with 
the inter-sensor distance set to 1 m and the local propa¬ 
gation speed of the wave equal to 1 m/s. The propaga¬ 
tion delay ti to estimate is thus located at 1 s. A set of 

1.1. d. plane wave sources impinges on the sensor set. The 
auto-correlation function of a source is the sine function 
with a large bandwidth. As presented in Equation (A.l), 
each correlation is normalized by its power. Each correla¬ 
tion is then degraded by some additive white sequence of 
Gaussian noise with variance set to 0.01. In Fig. □-Eni 
the top quadrant displays the estimated Green’s correla¬ 
tion. The distribution of the incidence angles that leads to 
this estimation of the Green’s function is displayed in the 
bottom-left. The color bar indicates the intensity map for 
the values of the distribution po. The estimated Green’s 
function is displayed on the bottom right. 

Note at this point that an alternative estimator for the 
propagation delay can be formed by retaining the largest 


7g(u,x,y) 


__ r 


l/ 27 r 



Figure 7: Setup (1). Uniform source field scenario. 

Number of sources = 10^ and bandwidth L. 
(bottom left) Incidence angles are uniformly 
distributed over the whole phase space, (bot¬ 
tom right) Provided that NL ^ 2ti, the differ¬ 
entiation of the Green’s correlation leads to the 
presence of a peak at 1 s. This is well matched 


with model (13) for a boundary free medium. 


7g(u,x,y) 



pe ~ ^^([0,27r]) 

Figure 8: Setup (2). Small number of sources scenario. 

Number of sources: 10^ and bandwidth L. 
(bottom left) Incidence angles are uniformly 
distributed over the whole phase space, (top) 
When NL < 2ti, the Green’s correlation is 
not properly estimated, (bottom right) As a 
consequence, the differentiation of the Green’s 
correlation is very inconsistent with model 
and thus denies propagation delay estimation. 


time difference of arrivals. This estimator is only valid in 
this setup which is the reason why it is not studied here. 

From the conducted simulations, uniformity of the dis¬ 
tribution of incidence angles coupled to a large number of 
sources appears essential to achieve passive propagation 
delay estimation. A slight deviation from the uniform dis¬ 
tribution leads to biased estimations. 

4.2.2. Performances study 

Two observation models were derived for the purpose 
of propagation delay estimation. The model in Equation 
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Figure 9: Setup (3). Band-limited source field scenario. 

Number of sources: 10^ and bandwidth L. 
(bottom left) As a special case from scenario 
in Fig. incidence angles are uniformly dis¬ 
tributed in a small subset of the phase space. 
This kind of setup is similar to that observed 
by m and m in a basin in California where 
the impinging directions of waves were con¬ 
fined to a subset of the whole phase space, 
(bottom right) Here, a bias B{ti) = t \/2 is 
created in the propagation delay estimation. 
Note however that if ^ = 0 is contained in the 
distribution, the estimation of the propagation 
delay in this setup may still be achievable. 
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0 time (s) 
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PB ~ n([0, 2 -k]) 



time (s) 


Figure 10: Setup (4). Non uniform source field sce¬ 
nario. Number of sources: 10^ and band¬ 
width L. (bottom left) Incidence angles 
are non uniformly distributed over the whole 
phase space, (top) The function that is es¬ 
timated in place of the Green’s correlation 
can be seen as a space-filtered version of it. 
Thus, the estimated Green’s function pre¬ 
vents parameters extraction from the applica¬ 
tion of estimators introduced for time-filtered 


(19) is specific to the passive approach as it focuses on the 
Green’s correlation function. From these models, two the¬ 
oretical lower bounds for the variance of the propagation 


delay estimation were derived, in Equations (20) and (16). 


The passive propagation delay estimation performances 
are computed based on setup ( 1 ), with a source field that 
has a varying bandwidth. The mean standard deviation 
for the passive and the active estimation of the propagation 
delay are compared to the theoretical bounds, see Fig. [m 



Figure 11: Standard deviation for the propagation de¬ 
lay estimation, (thick lines) (a) Performances 
for the passive propagation delay estimation 
(a’) Theoretical lower bound in the passive 
case, from the model of the Green’s correla¬ 
tion, in Equation ( [2Q| ). (dotted lines) (b) Per¬ 
formances obtained in the active case, based 


on model (15). (b’) Theoretical lower bound 
from Equation ( p^ , applicable for both the 
passive and the active case. 

The theoretical CRLB in the passive context is 2 decades 
lower than the actual standard deviation that was obtained 
in simulation. This discrepancy was expected, since the 


bound (20) was obtained from the observation model (19), 


which describes the Green’s correlation function. The 
propagation delay was not directly estimated from this 
data. Instead, the Green’s function was first retrieved by 


applying Equation (10) and only then was the propagation 


delay extracted. Thus, as quantified in [Appendix A[ huge 
performances losses were expected. The performances in 
the active context are sensibly better than in the passive 
context but both reach an experimental lower bound as in 


sources only, such as the one in Equation (f^). real-data experiment presented in section |T3 


(15) is a band-limited pulse in some noise, which is ap¬ 
plicable to the estimated Green’s function — whether in 
an active or a passive context. The model in Equation 


5. Perspectives and conclusions 

In this article, identification of a linear wave propaga¬ 
tion medium was introduced in a passive context, i.e. us- 


8 




































ing sources of opportunity only. The ward identity was 
extended to account for the spectral properties of acoustic 
sources. Estimators for the Green’s correlation function, 
the Green’s function and an example of parameter, the 
propagation delay, were then studied in the light of the 
temporal statistics of the source field. Performances from 
real-data experiments and simulations were compared to 
the theoretical bounds introduced in this paper. 

In the studied experimental setup, detection errors for 
the passive propagation delay estimation were shown to 
be less likely when the acoustic signals are lowpass filtered 
below 7kHz to lOkHz. In this case, performances were 
bounded by the temporal resolution due to the finite sam¬ 
pling rate. The accuracy of the inter-sensor distance esti¬ 
mation is comparable to the size of the sensors (~ 1 cm). 

Based on a geometrical analogy, an observation model 
with plane waves was derived for the passive context, al¬ 
lowing to account for the spatial properties of the source. 
It was shown that non uniformity of the propagation di¬ 
rections of the field, over the whole phase space, generally 
leads to severe estimation errors. This analogy leads to 
an observation model on the Green’s correlation function 
and the derivation of the accompanying Cramer-Rao lower 
bound. The passive and active estimation protocol were 
found to have comparable performances. 

Finally, it seems very beneficial to develop a technique 
that would perform passive propagation delay estimation 
from the Green’s correlation function, since the theoretical 
bound is attractively lower than the obtained experimental 
results. Furthermore, a model-based approach would allow 
to estimate parameters with a resolution that would rely 
less on the sampling rate, which was the limiting factor for 
the performances in the studied experiments. 

An extension with great potential in Structural Health 
Monitoring includes the passive localization of embedded 
or burrowed sensors and starts by extending the Ward 
identity to the concerned medium {e.g. concrete). 
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Appendix A. Passive identification protocol 

Measurements at locations x and y are times series x[n] 
and y[n]. Acquisition is assumed sampled at frequency 
1 /T Hz. Quantities of interest are the pressure ampli¬ 
tudes due to some sources of opportunity, that are seen by 
both sensors. Measurement noises are assumed to be i.i.d. 
zero-mean Gaussian random sequences, Wa,(t) and w^(t) 
respectively. The observation model is thus 

x[n] = /(nr,x) + Wa;[n] 

y[n] = f{nT,y)+vfy[n] 

In practice, the Green’s correlation 7 ^(t) in Equation <§ 
is estimated from finite-length and noisy signals. Expec¬ 
tations over time and space can be computed by splitting 
the signals x[n] and y[n] into K segments using a sliding 
window of constant duration. The average of all windows 
allows to estimate the spatial expectation, seen as an en¬ 
semble average. The kth window is thus the measurement 
between the starting index and the ending index 
where \bk — ak\ is constant for all k. Note that overlap be¬ 
tween windows is possible. An estimation of the Green’s 
correlation, can be 


x[n]^^^ _ — x[n — A simple proof leads 

to the conclusion that the variance of the ath order dif¬ 
ferentiation of x[n] is where (^) is the binomial 

coefficient. The order of the differentiation is given by the 
dissipation model. In the case of weakly viscous damping, 
(a = 3, see Equation Consequently, the variance of the 
estimation residue of the Green’s function is 20 times the 
variance of the estimation residue in the Green’s correla¬ 
tion. In comparison, the constant damping model used in 
geosciences leads to a factor 2 instead. As a consequence, 
although differentiation is a necessary step in the passive 
estimation protocol, if one wants to retrieve the Green’s 
function, one has to bear in mind that the variance of the 
estimation error will always increase, and that it depends 
on the damping model. 

Appendix B. A geometrical analogy 

The relationship between the time difference of arrival 
and the inter-sensor propagation is a function of the source 
position, see Eig. In a passive context, the position of 
the source is uncontrolled. We would like to link the time 
difference of arrival A that is induced by any source to an 
equivalent incidence angle Oeq of a plane wave source that 
would preserve exactly the same time difference of arrival 

A(s) = A(6»eq) = h COsiOeq), 

The position of the source is allowed to change while 
still preserving the time difference of arrival 

Vs, I d{s, x) — d{s, y) \ = vA. 




bk-n 


^ ^ i^k - Cik)V^ ^ 


^[l]y[l - n]. (A.l) 


In Equation (A.l), the correlation estimator is biased 


but has low variance m- The bias is convolutive (Bartlett 
window) but does not impact the position of the max. 
The emphasis is then put on the variance which has to 
be kept as low as possible. The normalizing constant 
Ek = power associated to the 

fields in the kth window. This normalization is important 
since inequalities on the energies of the fields can be due 
to inequal power of sources, geometrical attenuation and 
dissipation, but do not necessarily imply uncorrelatedness 
of the fields. 

By the application of the Ward identity (|^, the Green’s 
function is retrieved by differentiating the correlation ob¬ 


tained in Equation (A.l). Differentiation can be numeri¬ 


cally computed using different techniques. Whatever the 
technique, it is known that the variance of signal differ¬ 
entiated greatly increases. To provide an order of magni¬ 
tude, consider the basic case of differentiating using sub¬ 
tractions. Let x[n] be some zero-mean i.i.d. Gaussian ran¬ 
dom sequence. Its variance is The a-th differentiation 
of x[n] can be numerically computed using the formula 


where v is the wave velocity and (i(a, b) is the Euclidean 
distance between locations a and b. This relationship natu¬ 
rally casts the set of the positions allowed into a hyperbola 
1-L with the reduced equation 

The hyperbolas intersect the [x; y] segment at a = A/2, 
which then allows to find b = ^y^'(ut7^rAy(ut7^^~^^ 
Eventually, recall that hyperbolas have asymptotes with 
slopes a = ±b/a^ which defines the incidence angle of the 
equivalent plane-wave source 

6 eg = arctan {^\/{vh + A){vti - A)). 

As a conclusion, any time difference of arrival can be 
infered by a unique equivalent plane-wave source. 
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